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Abstract -We study numerically the wavefunction evolution of a Bose-Einstein condensate in a 
Bunimovich stadium billiard being governed by the Gross-Pitaevskii equation. We show that for 
a moderate nonlinearity, above a certain threshold, there is emergence of dynamical thermaliza¬ 
tion which leads to the Bose-Einstein probability distribution over the linear eigenmodes of the 
stadium. This distribution is drastically different from the energy equipartition over oscillator de¬ 
grees of freedom which would lead to the ultra-violet catastrophe. We argue that this interesting 
phenomenon can be studied in cold atom experiments. 


Introduction. — Starting from the famous Fermi- 
Pasta-Ulam problem m the interest to understanding 
of thermalization process in dynamical systems with a fi¬ 
nite number of degrees of freedom is continuously growing 
(see e.g. [3j[4]). At present the experiments with cold Bose 
gas in atom traps and optical lattices open access to ex¬ 
perimental investigations (see e.g. [5j 6j) stimulating the 
theoretical and experimental interest to this phenomenon 
(see e.g. (Tpl). 


The numerical analysis of dynamical thermalization 
in disordered nonlinear chains has been started recently 
showing that the quantum Gibbs distribution can appear 
at a moderate nonlinearity contrary to usually expected 
energy equipartition over linear modes (9 10 . Thus it is 
interesting to analyze the dynamics of a Bose-Einstein con¬ 
densate (BEC), described by the Gross-Pitaevskii Equa¬ 
tion (GPE), in a chaotic billiard where the quantum evolu¬ 
tion corresponds to a regime of quantum chaos and energy 
levels statistics described by the Random Matrix The¬ 
ory 11 13]. As an example we use a de-symmetrized 
Bunimovich stadium where the classical dynamics is fully 
chaotic (see e.g 


Model description. — The model is described by the 
GPE for BEC in the de-symmetrized Bunimovich stadium 
billiard with Dirichlet boundary conditions: 

= - 2- A^(f, t) + ip(r, t)\ 2 i>(r , t) (1) 

where we consider Ti = 1, m = 0.5. The height of the 
stadium is taken as h = 1 and its maximal length is 
l = 2 (see Fig. 0 . Thus the average level spacing is 
A « 47r/A « 7.04 where A is the billiard area. At /3 = 0 
the numerical methods of quantum chaos allows to de¬ 
termine efficiently about million of eigenenergies of linear 
modes and related eigenmodes 17 . For comparison, we 


12 


14 


also consider the case on a rectangular billiard with ap¬ 
proximately the same area as for stadium and with the 
golden mean ratio l/h = (1 + \/5)/2, h = 1. We note that 
for model 0 the spectrum of Bogoliubov excitations of 
BEC in a Bunimovich stadium had been studied in 18], 
but the question of dynamical thermalization has been 
not addressed there. We also point that the model 0 is 
described by the partial differential equation (continuous 
variables) thus being significantly more complex than the 
case of nonlinear chains studied in |9,10 


Indeed, 


We note that the chaotic opti- the prove of the existence of solution in ([T]) is a nontrivial 


cal billiards, created by a laser beam and containing cold 
atoms, have been already studied experimentally 15 16 


and hence our model can be investigated experimentally. 


mathematical problem which still remains open (e.g. the 
ultra-violet catastrophe would imply the absence of solu¬ 
tion). In the following we restrict our analysis to the GPE 
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case not going beyond this description. 


m’=5 m’=10 



Fig. 1: (Color on-line) Time evolution of wavefunction ampli¬ 
tude y , t) | in coordinate representation for an initial state 
of linear eigenstate mode m — 5 (left column) and m! — 10 
(right column) shown at t = 0 (top panels). Middle and bot¬ 
tom panels show the snapshots of corresponding distributions 
at t = 4 and t — 40 respectively. Here (3 — 10, color bars give 
values of \ip(x,y,t)\. 


The time evolution of 0 is integrated by small time 
steps of Trotter decomposition of linear and nonlinear 
terms with a step size going down to At = 4 x 10 -5 
to suppress nonlinear growth of high modes. We use 
N c = 1076(1085) linear modes of stadium (rectangu¬ 
lar) for linear part of time propagation doing the non¬ 
linear step with (3 term in the coordinate space with 
N p = 11207(12816) points inside the billiard. The lat¬ 
tice points are given by 79 x 79 = 6241 equidistant x-y 
coordinates for the square part of the stadium billiard 
(rectangular billiard), and rays with equidistance in an¬ 
gles for the circular part. The change of basis from co¬ 
ordinates to energies (and viceversa) is given by a uni¬ 
tary matrix in double precision. Similar to 19 , a special 


aliasing procedure is used with an efficient suppression of 
nonlinear numerical instability at high modes. This in¬ 
tegration scheme exactly conserves the probability norm 
providing the total energy conservation with an accu¬ 
racy better than 2% at largest value /3 = 20 and better 
than 1% at lower /? values. At any step the wavefunc¬ 
tion is expanded in the basis of linear modes </> m so that 
ip(x,y,t) = J2 m Cm(t)(l>m(,x,y). The averaging of proba¬ 
bilities w m (t) = \C m (t)\ 2 {Yj m w m = 1) over time gives 
the average probability distribution p m = (\C m \ 2 )t- 

We note that the quantum evolution of GPE 0 has 
been studied in the frame of quantum turbulence for a 
rectangular billiard 20 and for a 3D-cube [21]. However, 
in these studies there is energy injection/absorption at 
low/high modes to generate the Kolmogorov energy flow 
in space of linear modes (see e.g. 22 ,[23]). We also note 
that the time evolution of wave packet for the GPE in 
Bunimovich stadium had been simulated in [24] but only a 
spacial distribution had been considered there. In contrast 


to these studies we consider only unitary or Hamiltonian 
evolution given by 0 being interested in the distribution 
properties of probabilities p m over linear eigenmodes. In 
this respect our approach is different from other studies 
where the analysis had been concentrated on space fluc¬ 
tuations (see e.g. 25]). Also, as we will see below, there 
is a significant difference for the GPE evolution in chaotic 
and rectangular billiards. 

Time evolution. — Examples of time evolution for 
two initial eigenmodes m! = 5,10 are shown in Fig. [l] 
at j3 = 10 (video is available at 26]). They show that, 
due to nonlinearity inside the stadium, there are complex 
irregular oscillations of wavefunction with time. 

Another representation is obtained by considering the 
time evolution of probabilities w m (t ) in the basis of linear 
modes shown in Fig. [2] At moderate value (3 = 2 the 
probability remains located only in a few modes without 
thermalization and spreading over many modes. For larger 
value of /? = 10 the nonlinear spearing over modes goes 
in a more efficient way with many excited modes. Thus 
the dynamical thermalization is expected to be absent at 
small or moderate (3 < (3 C ~ 1 while for large nonlinearities 
/3 ~ 10 > /3 C we may expect the emergence of dynamical 
thermalization. 

For the initial state m! — 10 we have the approxi¬ 
mate energy value E rn « m'A & mv 2 / 2 70, where 

v ~ y/280 is a velocity of classical particle and a time inter¬ 
val between collisions is approximately r co i ~ h/v ~ 0.06. 
Thus during the time t = 40 we have approximately 
N co i = t/r co i ~ 670. Dynamical thermalization is reason¬ 
ably achieved for time interval t G [20,40] as it is visible 
in middle and bottom panels of Fig[2] where w m (t) have 
thermalized like distribution with f3 = 10. 

Bose-Einstein thermal distribution. — To char¬ 
acterize the dynamical thermalization in more detail we 
assume that a moderate nonlinearity acts as a certain ther- 
malizer which drives the system to a thermal equilibrium 
over quantum levels of the stadium. At the same time we 
assume that the nonlinear term is not very strong so that 
it does not affect significantly the average linear eigenen- 
ergies. Indeed, on average we have P\^\ 2 ~ /3/A « A 
for (3 ~ 10, so that indeed, the nonlinear energy shift is 
moderate at such values of /3. 

Thus we expect that nonlinearity generates a dynamical 
thermalization over the quantum billiard energy levels. In 
such a case we should have the standard Bose-Einstein 
distribution ansatz over energy levels E m [27 : 

Pm = l/[exp [(E m -E g - p)/T] - 1] (2) 

where E g = 13.25 is the energy of the ground state, T 
is the temperature of the system, p(T) is the chemical 
potential dependent on temperature. The values E m are 
the eigenenergies of the stadium at (3 = 0. The param¬ 
eters T and p are determined by the norm conservation 
Ylm=i Pm = 1 (we have only one particle in the system) 
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Fig. 3: (Color on-line) Entropy dependence on energy S(E) 
obtained from the GPE time evolution of initial linear eigen¬ 
states with 1 < m < 50 for the stadium (circles) and rectan¬ 
gular (crosses) billiards for nonlinearity /3 = 2, 5,10, 20 marked 
on each panel. Here the average is done over time intervals 
t G [4,5] (blue), t G [16,20] (green) and t G [20,40] (black). 
The red curve represents the Bose-Einstein ansatz © while 
the orange dashed curve shows the case of energy equiparti- 
tion over first 50 modes of the stadium. 


Fig. 2: (Color on-line) Time evolution of probabilities ic m (t) 
in linear mode basis for initial state m! — 10 at = 2 (top 
panel), m — 10,20 at (3 = 10 (middle and bottom panels 
respectively). The probabilities Wm(t) are averaged over time 
interval St = 0.2 to reduce fluctuations; color bar shows values 

of logio Pm(ra'). 


and the initial energy E m Pm = E. The entropy S 
of the system is determined by the usual relation 27 


S = - Yjm Pm In Pm- The relation © with normalization 
condition determines the implicit dependencies on tem¬ 
perature E(T), S(T ), /i(T). 

The advantage of energy E and entropy S is that both 
are extensive variables, thus they are self-averaging and 
due to that they have reduced fluctuations. Due to this 
feature S and E are especially convenient for verification 
of the thermalization ansatz. To check this ansatz we 
start from an initial linear mode m' which corresponds 
to the system energy E ~ E m t and follow the GPE time 
evolution of probabilities w m (t) determining the value of 
entropy S from obtained average probabilities p m . Con¬ 
sidering the initial states with 1 < m' < 50 we obtain the 
numerical dependence S(E) shown by symbols in Fig. [3j 
This dependence is compared with the analytic curve fol¬ 
lowing from the Bose-Einstein ansatz © which gives the 
dependencies E(T) and S(T) and hence provides the an¬ 
alytic dependence S(E) shown by the red curve in Fig. [ 3 ] 

The data of Fig. [3] show that even at large /3 = 20 there 


is no thermalization for the rectangular billiard. We at¬ 
tribute this to the fact that the ray dynamics is integrable 
in this billiard and thus it is much more difficult to reach 
onset of chaos for the GPE in this billiard at moderate 
nonlinearity studied here. 

The situation is different for the stadium: at f3 = 2 only 
a few modes are populated, at /3 = 5 the number of modes 
is increased but still the numerical data for S are very far 
from the thermalization red curve (at least on the time 
scale reached in our numerical simulations). However, for 
/3 = 10, 20 we find that the numerical data at large times 
t > 15 follow the theoretical curve S(E) given by the Bose- 
Einstein thermalization. A small visible deviation from 
theory is still visible since the numerical points are sys¬ 
tematically slightly below the theory curve. We attribute 
this to a finite computation time which apparently is not 
long enough to visit all regions of multi-configurational 
space with sufficiently large statistics. On the basis of 
obtained data we can conclude that the dynamical ther¬ 
malization in the stadium sets in for /? > /3 C « 7 « A. 
We also checked that the initial states, which represent 
a linear combination of a few eigenmodes, also follow the 
theoretical red curve in Fig. [ 3 ] at j3 > /3 C (e.g. two modes 
m = 10,15 at P = 10). 

Another way to check the thermalization predictions 
is to determine T from the numerical values of F, S 
which, according to the Bose-Einstein ansatz, give in¬ 
dependent values T\(E) and 72(5). The average value 
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Fig. 4: (Color on-line) Theoretical dependencies given by Bose- 
Einstein ansatz j2| and shown by the red curves for T(E) (top 
left panel), T(S) (top right panel), p(E) (bottom left panel), 
p{T) (bottom right panel); the numerical results, obtained 
from GPE in the stadium, are shown by black points at = 10 
with averaging over the time interval [20,40]. For representa¬ 
tion convenience we show these dependencies using rescaling 
to maximal values of variables corresponding to initial state 
with ml — 50: E max — 414, Smax — 4.8, T max — 387.75, 

| Pmax | = 1500.45. 


T = (Ti +T2)/2 is shown in Fig.[4]as a function of numer¬ 
ical values of E and S. We see that the numerical points 
are in a good agreement with the analytic curves (apart 
of small systematic displacement of numerical points dis¬ 
cussed in the paragraph above). In a similar way we make 
a comparison between the theory and numerical data for 
dependencies p(E) and /i(T) shown in bottom panels of 
Fig.@ Again we find a good agreement between the nu¬ 
merical data and the Bose-Einstein thermalization distri¬ 
bution (J2|. 


The validity of the Bose-Einstein thermalization ansatz 
© leads to a striking paradox pointed already for nonlin¬ 


ear chains in 10 : formally the GPE in stadium gives a 


system of equations for nonlinear coupled oscillators (we 
have nonlinear coupling between linear oscillator modes of 
the stadium) with a moderate nonlinearity. The usual ex¬ 
pectations of the statistical mechanics predict the energy 
equipartition between these modes [10,27 . If the number 
of modes is infinite then we should have ultra-violet catas¬ 
trophe with probabilities p m ~ 1/E m at high modes m and 
the global temperature approaching zero as T ^ 1 /m m a X 
for initial excitation with m' ~ 10 (here m max is the max¬ 
imal mode number). This classical thermalization ansatz 
for mmax — 50 is shown by a dashed curve in Fig. [3] and 
the data clearly show that it is very different from the nu¬ 
merical data which are close to the Bose-Einstein ansatz 
(we note that numerically the quantum Gibbs distribu¬ 
tion 110] over quantum levels of stadium gives the results 




0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 

Fig. 5: (Color on-line) Time averaged probabilities pm(jn) at 
stadium eigenstate m for initial state ml \ the time averaging 
is done for time intervals [20,40]; the panels show data for 
1 < m',m < 30 in x,y axes respectively. Here we show the 
cases: f3 = 2 (top left panel), [3 = 5 (top right panel), [3 = 10 
(bottom left panel), the theoretical Bose-Einstein distribution 
|2]) (bottom right panel). The values of pm(jn) are shown by 
color with the corresponding color bars for each panel. 


being rather close to those of © since at low temper¬ 
atures and large m values both distributions are rather 
similar). Thus our data clearly show the emergence of the 
dynamical thermalization described by the Bose-Einstein 
distribution in a chaotic billiard for moderate nonlinearity 
(3 > /3 C ~ A. 

A more detailed check of the Bose-Einstein distribu¬ 
tion requires a direct comparison of numerically obtained 
probabilities p m (ra') with the theoretical expression © 
for each initially excited mode ml. We show such data 
in Fig. [5] for [3 = 2, 5,10. It is clear that there is no 
thermalization at (3 = 2,5 since a large fraction of prob¬ 
ability remains at the initially populated state m = ml . 
For [3 = 10 we see that the probability at initial state 
ml drops significantly indicating emergence of dynamical 
thermalization. However, still the numerical probabilities 
at large m have larger values compared to those of the 
theory Q shown in the bottom right panel of Fig. [ 5 ] We 
attribute this to the fact that our total computation time 
is not large enough to have good statistical data for av¬ 
erage values of p m (ra') which require good averaging and 
long computation times. Such a problem had been visible 
in the numerical simulations with nonlinear chains [9 
where the time of simulations have been by a few orders 


10 


P-4 






































Dynamical thermalization of Bose-Einstein condensate in Bunimovich stadium 


of magnitude larger than here. At the same time the ex¬ 
tensive property of energy E and entropy S makes them 
self-averaging and more stable in respect to fluctuations 
thus allowing to compare them with the theory 0 at sig¬ 
nificantly shorter time scales. 

Unfortunately the large scale simulations of the GPE for 
stadium are rather heavy and time consuming since they 
require transformations from coordinate to linear mode 
space and small integration time step with aliasing proce¬ 
dure to suppress numerical instabilities of high modes. It 
is possible that the numerical codes can be improved al¬ 
lowing to reach larger time scales but this requires further 
studies going beyond the scope of this work. 

Finally we discuss a preparation of one or a few ini¬ 
tial eigenstates considered above for the time evolution 
and dynamical thermalization. It is clear that the ground 
state of the billiard is relatively easy to prepare since it 
is the final state in a process of relaxation and also since 
it is compact in space being close to a coherent state of 
a harmonic trap. An excited state can be produced from 
the ground state applying a monochromatic driving (oscil¬ 
lation) of the billiard that creates an effective ac-potential 
V ac = fxcos(u:t ) if one goes to the oscillating frame (see 
e.g. [28]). In a chaotic billiard dipole matrix elements have 
transitions between all energy eigenstates 28 and thus a 
resonant transition will populate one or a few states being 
close to the resonance E n « Eq + Two. We note that such a 
method demonstrated already its efficiency for excitation 
of high energy states for chaotic Rydberg atoms (see e.g. 


29 30 ). 


standing. We hope that our results will stimulate further 
research in this field of fundamental aspects of nonlinear 
dynamics and thermalization onset in systems with large 
but finite number of degrees of freedom. 

The modern progress in the cold atom experiments al¬ 
lows to investigate a dynamics of Bose gas and Bose- 
Einstein condensates [5] 6] while the chaotic billiard for 
such atoms can be created by optical beams [ISpC . Thus 
we think that the model 0 can be realized with cold atom 
experiments. 





Discussion. — Our studies of the GPE in the Buni¬ 
movich stadium billiard show that for a moderate nonlin- 
eariy parameter above a certain threshold /3 > /3 C « A 
the nonlinear Hamiltonian dynamics leads to emergence 
of dynamical thermalization over the linear billiard modes 
which is well described by the Bose-Einstein distribution. 
This distribution is strikingly different from the usually 
expected energy equipartition over modes 


10 27 which 


would lead to a violet catastrophe with a significant prob¬ 
ability transfer to higher and higher modes of the chaotic 
billiard. The established validity of the Bose-Einstein dis¬ 
tribution, together with the previous studies of dynami¬ 
cal thermalization in nonlinear chains (9 10 , leads to an 
unexpected conclusion about emergence of quantum dis¬ 
tributions over linear energy modes in systems of coupled 
nonlinear oscillators at moderate nonlinearity. This result 
is drastically different from the standard energy equiparti¬ 
tion picture expected for nonlinear dynamics of oscillator 
systems 1-4 27 . 


The described picture of “quantum” dynamical ther¬ 
malization for the GPE in chaotic billiards requires a bet¬ 
ter understanding of nonlinear dynamics in systems with 
many degrees of freedom. It is known that slow chaos, like 
the Arnold diffusion 31-33 , and an anomalous diffusion in 


disordered nonlinear chains (see e.g. 19,34-36]) generate 
a number of features which still wait their deep under¬ 



x x 


Fig. 6: (Color on-line) Poincare sections (x,p x ) at y — 0.5, p y > 
0. Left and right columns correspond to energy E — 2 and 
E — 18 respectively. Top panels show the entire Poincare 
sections (10 chaotic orbit up to time t < 10 4 ) and middle panels 
show zoom marked in top panels (adding 15 trajectories in the 
integrable region to time t < 10 4 ), one invariant curve of each 
panel is highlighted with red (gray) color inside stability island 
at middle panels. Bottom panels show dynamics in (x, y) plane 
with stable orbits from middle panels (same red color) and 
chaotic orbits (orange color also shown in top panels) up to 
times t — 10 2 ; dashed horizontal lines mark y — 0.5 used for 
the Poincare sections. 

Another promising possibility can be an experimental 
realization of a harmonic Sinai billiard, or Sinai oscillator. 
An example of such a billiard is described by a classical 
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Hamiltonian H = ( p x 2 + x 2 )/2 + ( p y 2 + 2y 2 )/2 with a 
rigid disk of radius ry — 1 located at x == y = — 1 (thus 
the ratio of frequencies in x and y is irrational). The har¬ 
monic potential can be realized by optical traps while the 
repulsive rigid disk can be created inside by an additional 
laser beam with such a frequency detuning that it acts as 
a repulsive potential for cold atoms. Examples of typi¬ 
cal trajectories in such a billiard are shown in Fig. [6] In 
the same figure we also show the Poincare sections con¬ 
structed in (x,p x ) plane at y = 0.5, p y > 0 with energies 
E — 2 and E m 18, when the size of oscillations of atom 
is larger than r^. In this regime almost all phase space 
is chaotic (the domains of integrable dynamics are very 
small). Thus we think that the GPE in such a harmonic 
Sinai billiard will show all the effects of dynamical ther- 
malization discussed above for a more convential case of 
the Bunimovich stadium. We expect that such a system 
can be more simple for experimental investigations. Also 
in such a billiard a coherent state of the harmonic potential 
can be created experimentally and can be used as an ini¬ 
tial state with energy being close to the energies of linear 
eigenmodes of such a billiard. We expect that experimen¬ 
tal investigations of the GPE in a harmonic Sinai billiard 
(or in the Bunimovich stadium) will allow to understand 
the fundamental aspects of dynamical thermalization. 

We thank D.Guery-Odelin for useful discussions of cold 
atom physics. 
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